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ABSTRACT 

We investigate fluid motions near tlie midplane of vertically stratified ac- 
cretion disks with highly resistive midplanes. In such disks, the magnetorota- 
tional instability drives turbulence in thin layers surrounding a resistive, stable 
dead zone. The turbulent layers in turn drive motions in the dead zone. We 
examine the properties of these motions using three-dimensional, stratified, lo- 
cal, shearing-box, non-ideal, magnetohydrodynamical simulations. Although the 
turbulence in the active zones provides a source of vorticity to the midplane, no 
evidence for coherent vortices is found in our simulations. It appears that this 
is because of strong vertical oscillations in the dead zone. By analyzing time 
series of azimuthally- averaged flow quantities, we identify an axisymmetric wave 
mode particular to models with dead zones. This mode is reduced in amplitude, 
but not suppressed entirely, by changing the equation of state from isothermal to 
ideal. These waves are too low-frequency to affect sedimentation of dust to the 
midplane, but may have significance for the gravitational stability of the resulting 
midplane dust layers. 

Subject headings: protoplanetary disk,magnetohydrodynamics 
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Introduction 



A long standing problem in star formation concerns the accretion of high-angular mo- 
mentum material in disks around protostars. The rediscovery of the magnetorotational 
instability (MRI) in an astrophysical context, along with the discovery of the favorable 



transport properties of the turbulence that results from it (Hawley et al. 1995) led to the 



widespread belief that a solution to this puzzle had been found. However, the instability 
requires a critical coupling strength between field and fluid ( Blaes fc Balbus||1994 Jin||1996 ). 
Noting the high densities and low temperatures likely at the midplanes of protoplanetary 
disks, Gammie (1996) developed a model including a weakly-coupled, MRI-stable "dead 



zone" at the midplane of such disks. Since this prediction, considerable effort has gone into 
elucidating the conditions under which such a model might exist. Most of this work has 
centered on chemical studies of the conditions under which such a dead zone might form 
and what its radial and vertical extent might be. Recently, a number of numerical studies. 



pioneered by Fleming & Stone (2003), have begun to explore the dynamical consequences 



of the dead zone model Fromang & Papaloizou (2006). Turner et al. ( |2007 ) have developed 
an interesting coupled chemical and magnetohydrodynamical (MHD) model in which the 
dead zone disappears because of vertical mixing of ionized gas to the midplane. In their 
picture, the dead zone is absent in the sense that large scale magnetic fields with positive 
Maxwell stress, {—B^B^), form in the midplane and transport significant amounts of angular 
momentum even though the midplane is still MRI stable. 

Aside from acting as a conduit for material flowing onto the protostar via accretion, 
protoplanetary disk midplanes are also the sites of planet formation. One of the outstanding 
issues in planet formation is the collection of micron-sized dust into planetesimals, rocky 
objects of order 1 km in size. In this context, turbulent fluid motions more general than the 
correlated fluctuations that drive accretion become very important. This is because drag 



forces couple the dust and gas. Fully MRI-active simulations, Johansen & Klahr (2005) 



find preferential concentration of dust in turbulent over-pressures. However, these same 
turbulent motions increase the velocity dispersion of the dust grains, possibly leading to de- 



structive collisions. The dust-gas coupling can drive strong instability and clumping ( Youdin 



& Johansen 2007 Johansen & Youdin 2007), which, in combination with MRI turbulence. 



allows direct formation of planetesimals by gravitational instability (Johansen et al. 2007). 



Together, these results suggest that turbulence plays a complex role in planetesimal forma- 
tion. As a first step to considering more complex dust-gas models including dead zones, it 
is critical to characterize the fluid motions in the MRI-stable midplane. 

Here, we address two specific questions about the fluid motions in dead zones: do 
coherent, two-dimensional (2D) vortices form in dead zones? And, what are the wave-like 
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motions seen in our previously published models of dead zone dynamics (e.g. Oishi et al. 
2007)? Vortices have been studied both for their ability to drive purely hydrodynamic 



accretion 



Umurhan & Regev (2004) and their ability to trap dust particles and thus act as 



sites of enhanced planetesimal formation ( Bracco et al.||1999 Johansen et al.||2004 ). Johnson 



& Gammie (2005b, hereafter JG05) studied vortex dynamics in a 2D, compressible, shearing 
sheet model and concluded that random vorticity perturbations do form coherent large-scale 
vortices, but that energy in such structures decays away as They also suggested that 

turbulent overshoot from the active layers might provide a vorticity source at the midplane. 



Using three-dimensional (3D), compressible but unstratified simulations, Shen et al. (2006) 



showed that small amplitude vortices are unstable to elliptical instabilities (Kerswell 2002) 



and that finite-amplitude trailing waves are torn apart by a Kelvin-Helmholtz like instability 
before they have a chance to form coherent vortices. The anelastic simulations of |Barranco &: 



Marcus (2005) show large, coherent vortices forming at 1-3 H from the midplane, where H 



is the disk scale height. This unusual and quite unexpected result suggests that stratification 
plays an important role in disk vortex dynamics. Because our time-steps are Courant-limited 
by the Alfven velocity, we restrict our domain to \z\ < 2H to avoid low-density regions with 
extremely high Alfven speeds. Nonetheless, we study the effect of stratification on the flow 
inside and outside of the dead zone with the goal of understanding if MRI turbulence in 
active layers can sustain vortices. 

In section [2] we briefly describe our numerical techniques and how we tested them, while 
section [3] contains the results of our 3D simulations, and section |4] reviews aspects of vortex 
dynamics and presents relevant 2D simulations. Finally, we discuss our results and present 
conclusions in section [H 



2. Numerical Method 



2.1. Model Description 



We use the Pencil Code (Brandenburg & Dobler 2002) to solve the resistive MHD 



equations in the vertically stratified, shearing box limit (Hawley et al. 1995). 



All models are 3D, isothermal, non-ideal MHD, except for one ideal MHD control run. 



and one non-ideal, non-isothermal run (see section 3.3.1). They are run on a domain of 
IH X AH X 4if, where H is the disk scale height. The numerical method is stabilized with 
sixth-order hyperviscosity, hyperresistivity, and a hyperdiffusive operator on the density 



as described in Oishi et al. (2007). All the models are initialized with a MRI unstable 
zero-net flux initial magnetic field given by Bz{x) = Bosin{27ix / L^) with initial plasma 
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P = 2/io-Pg/-Bo — 400. The domain has periodic boundary conditions in y and z, and 
shearing periodic boundaries in x. 

The non-ideal models use one of two ^-dependent resistivity profiles. One is motivated 
by balancing cosmic ray-ionization to a critical depth in column density given by [Fleming fc 



Stone (2003). These models were previously described in Oishi et al. (2007) and are hereafter 



referred to as "FS runs" . These models were run with cubical zones of size 32 and 64 zones 
per scale height. Unless otherwise noted, we use a resolution of 64 zones per scale height 
for the FS runs. The second profile, described below, uses hyperbolic tangent functions and 
models using it were run with 32 zones per scale height. 



2.2. Hyperbolic Tangent Resistivity Profile 

This alternate profile offers a sharper transition between active and dead zones than the 
FS profile, and allows the thickness of the dead zone to be controlled independently from 
its depth. Though simpler, it is less physically motivated than the FS profile. The profile is 
given by 

where tjq is the midplane resistivity, zq is the transition height and Az is the width of the 
transition region, set in this work to Az = 5dx, where grid zones have size dx. We set 
rjo = 1.67 X 10^^, corresponding to magnetic Reynolds number Rcm = c^/iVmid^) = 30 at 
z < \zq\. This gives all our tanh models the same value of Rcm in the dead zone as our 
fiducial FS run has at z = 0. 

The Maxwell stresses as a function of height for both sets of profiles are given in Figure [T] 
The figure emphasizes the main utility of the tanh profile: because the dead zone resistivity 
is independent of its size, we can study larger dead zones at a given resolution. 



2.3. Tests 



In the numerical study of non-axisymmetric perturbations in disks, a series of test 
problems are becoming standard. The vortical and non-vortical shearing waves have closed 



analytic solutions (in the former case, a WKB approximation) derived by Johnson & Gammie 



(2005a). These waves have been used by JG05 to test a code derived from ZEUS (Stone & 



Norman||1992 ) and the higher-order Godunov code ATHENA ( Shen et al.||2006 ). Aside from 
providing an analytic test solution for comparison, the incompressible wave is particularly 
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useful for determining the amount of aliasing present in the numerical scheme. As a wave 
swings from leading to trailing, it wraps toward axisymmetry. However, in doing so the 
number of zones resolving the wave must drop. When the number of zones resolving each 
wavelength drops below some threshold, the code may unphysically transfer ( "alias" ) power 
from the trailing wave into a new leading wave that will again swing amplify. This danger 



is highlighted by both Umurhan & Regev (2004) and Shen et al. (2006), as such power may 



spuriously drive coherent vortices and angular momentum transport. 

The Pencil Code uses spatially-centered finite differences and Runga-Kutta timestep- 
ping, meaning that its discretization lacks any formal algorithmic viscosity. Thus, it should 
perform well on this test. Ultimately, however, this property works against the code: without 
an explicit viscous term, it does not dissipate energy that moves to smaller scales, where an 
unphysical buildup occurs. In this sense, the shearing waves mimic the effect of a turbulent 
cascade: as time progresses, kx{t) decreases and the waves become less and less resolved. 
This is not a turbulent cascade of course, but the code is nonetheless forced to resolve smaller 
and smaller structures, leading eventually to a crash. To avoid this in these tests, as in our 
actual runs, we use a sixth-order hyperviscosity scaled by uq for stability, setting the hyper- 
viscous Reynolds number at the grid scale Re^ = udx^/v^ ~ 1, where u is the maximum 
velocity on the grid. 

Figure |2] shows the incompressible shearing wave for runs with resolution = 64, 128, 256, 
demonstrating the ability of the Pencil code to accurately reproduce the analytic solution. 



The aliasing time in terms of the orbital frequency is (Shen et al. 2006) 



talias^O — — ( ; ) ; (2) 

q \ny ky J 

where n is an integer, q = dln^o/dlnr is the shear parameter, and Uy = kyLy/2'K is the 
dimensionless azimuthal wavenumber of the wave. The oscillations around t = 20 appear 
to be small amplitude sound waves present in our compressible solution but not in the 
incompressible analytic solution. 

In Figure |2| the = 64, 128, and 256 models run to tQo = 100, where Qo is the 
local orbital frequency. This is well beyond the expected aliasing time taUas^o — 45. 3n 
for Nx = 128 and 90. 6n for A^^ = 256. The figure makes an interesting point: aliasing 
does indeed occur for 64^, but because the hyperviscosity is a strongly non-linear function 
of resolution, the spurious energy injected by aliasing becomes trivial before the first 128^ 
aliasing time. Even at 64^, each successive aliased wave has less power, suggesting that 
sustained vortex activity is not likely to be powered by this numerical effect for long times. 
There is no aliasing at all at 128^ and 256^ resolutions, because the hyperviscosity damps the 
signal to machine precision before the first aliasing time. We have additionally performed a 
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convergence test using resolutions of 64, 128, and 256. Using hyperviscosity, the average error 
over tQ < 10 reaches a constant value at 128^ and does not further converge. However, when 
we stabilize the algorithm using a physical Laplacian viscosity instead of hyperviscosity, 
the error converges at roughly ~ where N is the resolution. We also 

consider the compressible shearing wave. Figure [3] demonstrates the code's ability to resolve 
the wave with small error to late times with moderate resolution. 

Finally, we successfully tested against the 3D, hydrodynamic, non-linear solution of the 



shearing box equations derived by Balbus & Hawley (2006). Figure shows the kinetic 



energy of the wave as function of time. In this test, as in the others, the wavenumber drops, 
so that at larger times the error necessarily grows. 

The incompressive shear wave test gives us a strict metric for aliasing: for the lowest 
resolution, we used H = 64dx, and saw trivial amounts of aliasing. Doubling the resolution 
to H = 128dx, we see no aliasing. Because the amount of vorticity increases in our MHD 
simulations when we increase the resolution from H = 32dx to GAdx (see figure [s]), we are 
confident that any sustained vorticity in our MHD simulations does not come from numerical 
aliasing, and that the code accurately tracks any vortical motions present. 



3. Results 

The excitation and saturation of the MRI occur by t/torb = 5 and the dead zone is 
clearly defined by orbit ~ 10 (see section 3.3 for details). Here, we will restrict our analysis 
to times later than 25 orbits in order to avoid transients. Figure [6] shows time series of 
vorticity (a; = V x u), and kinetic and magnetic energies. Although the saturated kinetic 
energy drops by over an order of magnitude between the Rcm = oo and the Rcm = 3 
models, the vorticity only drops by a factor of a few. This can be explained by the presence 
of residual vorticity in the dead zone. 



3.1. Vorticity and Flow Dimensionality 

The MRI produces significant amounts of vorticity, much of which is retained when a 
dead zone is introduced. Here, we investigate whether or not this vorticity can coalesce into 
coherent vortices that could be significant as a natural environment in which to trap dust 
and thus accelerate protoplanet formation. 
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In a thin accretion disk, the vertical gravitational acceleration is given by 



g 



(3) 



where Q is the orbital frequency. In the local approximation, we take Qq = Q{R) where 
R is the assumed central radius of the shearing box. For an isothermal disk in hydrostatic 
equilibrium, the Brunt- Vaisala frequency of buoyant vertical oscillation is 



gdp 

p dz 



niz^ 



(4) 



We can diagnose how stratified the flow is by comparing the vertical component of the 
vorticity ujz to the Brunt- Vaisala frequency in a version of the internal Froude number. 



e.g., Barranco &; Marcus 2005) 



Fr 



ujJ{2N). (5) 

vanishes at the midplane, and thus Fr formally diverges as 2; — > 0. Nonetheless, Fr provides 
a useful diagnostic for the degree of stratification in the flow: when Fr< 1, internal gravity 
waves rapidly homogenize vertical disturbances in the vortex, so the system is strongly 
stratified, and acts effectively as a 2D flow in which energy cascades to larger scales, similarly 
to the Earth's atmosphere. Thus, Fr is an effective measure of the dimensionality of a vortex 
in the disk plane. 

Previous results show rapid vortex growth in purely 2D systems integrated over z 
(Umurhan & Regev 2004 Johnson fc Gammie||200"5b ), and at Fr < 1 in 3D stratified 



systems (Barranco & Marcus 2005). Figure 7 shows the vertical behavior of Fr in our FS 



models. All but one of our simulations are purely isothermal, with P = cip, so gravity waves 
themselves are excluded. However, because the basic dynamical properties of the MRI are 
not very sensitive to the equation of state ( Stone et al.|1996 ), we do not expect that including 
gravity waves will change the Fr profile. We confirm this expectation in figure [7| which also 
includes our control Rcm = 30 run with a non- isothermal, 7 = 5/3 ideal gas equation of 



state (see section 3.3.1 for more details). The Fr profile is roughly similar to the run with 
isothermal Rcm = 30. 

In our ideal MHD run the flow is effectively 3D through almost the entire domain — the 
MRI is a strong enough vorticity source to overwhelm the homogenizing effect of vertical 
oscillations almost everywhere. On the other hand. Figure [7] suggests that the dead zone 
models have an effectively 2D flow at nearly all heights, most particularly in the dead zone 
of the Rcm = 3 model. Near the midplane, Fr diverges, as the stratification there goes to 
zero, so we would expect the very near midplane region to act as a 3D flow, consistent with 



the Barranco & Marcus (2005) results showing disappearance of their imposed vortices in 
that region. 
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However, Figure [8] shows that no vortices form in either the dead or the active zone. In 
the dead zone, the vorticity remains confined to elongated, nearly axisymmetric bands. In 
fact, no difference is seen in ojz between heights within the dead zone-only between dead and 



active zones. Thus, our results combined with Shen et al. (2006) and Barranco & Marcus 



(2005), show that even in the presence of a vorticity source such as active zones, stratification 



is necessary but insufficient to produce vortices. 

Table [T] shows the distribution of enstrophy among x,?/, and z components for FS models 



and the ideal MRI run. A strongly two-dimensional fiow shows a preferred direction (Bran- 



denburg et al. 1995). The ideal MRI run and the active zones of the FS models agree well 
with an isotropic distribution of vorticity. The dead zones show a strong preference for the y 
direction, indicating a circulation in the x — z plane, as indeed images bear out (see below). 
The coherent anticyclonic vortices we set out to investigate would have shown a strong z 
component. 

Figure [9] shows the components of velocity dispersion for each of the FS runs as a 
function of time, computed separately inside and outside of the dead zone. The z component 
of velocity dominates within the dead zone, while outside of it, it is a factor of a few smaller 
than the others. (Note that the oscillation in the dead zone velocities is real, but the plot 
was made with relatively crude time resolution, 5t = O.btorb and some aliasing from higher 
frequency signals is certainly present.) The dead zone motion is clearly not dominated by 
horizontal fiows. 



3.2. Morphology and velocity vectors 



Figure 10 shows images of density overlaid with velocity vectors for the midplane (the 



x-y plane at z = 0) and Figure 11 the x-z plane at ?/ = for each of the FS runs. In 
the dead zone midplanes (the right three panels), compressive waves dominate the velocity 
field: velocity vectors lie almost entirely orthogonal to the density striations in the x-y plane. 
However, in the x-z plane, inertial mode signatures appear: the velocity is dominated by the 
z component, traveling up and down in well defined vertical bands. 

The dead zone shows ordered alternating between positive and negative values across 



the dead zone that become more pronounced as Rcm decreases (Fig. 11, right three panels). 
This suggests that the vertical oscillations are at ^ = 0, with energy transport from gravity 
waves purely vertical between the active zones and the midplane. This is in stark contrast to 
the fully MRI-active case. Near its midplane, there are no well-defined vertical oscillations, 
and the Rcm = 100 model does not show nearly as well-ordered motion as the larger dead 
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zone models. This is perhaps not surprising given that the fastest growing MRI mode has a 
growth rate Qmri ~ ^O; which is greater than until around a scale height. 



3.3. Wave Modes 

The MRI is not operating in the densest parts of the disk in our dead zone models, so 
we expect that the motions excited in this region will take the form of linear wave modes 
stochastically excited by the active zone turbulence. The following analysis is motivated 
by the low-frequency oscillation in the volume-averaged kinetic energy most clearly shown 
in the Rcm = 3 data between 60 < t < 100 in Figure [6j Fourier analysis of this kinetic 
energy time series for runs with a dead zone show a clear peak in frequency space that shifts 
slightly to lower frequency with increasing dead zone size, suggesting the presence of a large 
amplitude coherent oscillation in these runs. 

In order to better understand these dead zone oscillations, we take temporal power spec- 
tra of the radial and vertical velocity components Ux, and u^, and the density perturbation 



5p = p{x,y, z) — po{z). Arras et al. (2006) and Brandenburg (2005) used similar methods to 



diagnose global modes of oscillation in fully magnetorotationally turbulent disks. The for- 
mer were able to isolate acoustic and inertial oscillations. However, the latter found no clear 
wave signatures for MRI turbulence, though they did recover acoustic and inertial modes for 
forced hydrodynamic turbulence in a shearing box. 



Figure 12 shows spectra computed by taking an FFT of a time series of each variable 
averaged in the y-direction at every (x, z) point on the grid of our low resolution model. The 
resulting power spectra were then averaged to raise signal-to-noise. Azimuthal averaging 



allows us to eliminate much of the MRI power in the active zones dArraset al.||2006D , though 



this may not be a significant source in some of our larger dead zone models. The dominant 
peak in all variables at rzi^ax — 0.23r2o shifts only very slightly, to lower frequency, as the 
dead zone is made thicker by a factor of two in the tanh runs. This suggests that the vertical 
thickness of the dead zone is not setting the oscillation frequency. We have also determined 
that this characteristic frequency is unaffected by the radial {x) extent of the box by running 
the Rcm = 30 model with = 2 and confirming that the oscillation frequency is unchanged. 

Because the box is periodic, standing oscillations can be excited, with discrete frequen- 
cies. The lowest frequency acoustic mode is in the vertical direction and has a frequency 
w ^ Cs2tt /Lz = l.llf^o, which is clearly larger than vumax- The isothermal equation of state 
used in these models precludes gravity waves (which would also have zu < Qq), and so we 
tentatively identify them as an inertial oscillation. Furthermore, the small shift seen in the 
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tanh runs as the dead zone thickness increases is actually to lower frequency, while a higher 
boundary would drive internal oscillations at a higher frequency (see Eq. |4]) . 



3.3.1. Non-isothermal models and Buoyancy Forces 

The low frequency of these oscillations is suggestive of a gravity mode. However, these 
runs are strictly isothermal, P = c'^p with Cs constant everywhere, precluding buoyant re- 
sponses and thus gravity waves. 



Recently, Bai & Goodman (2009) noted that isothermal equations of state might spu- 
riously enhance vertical mixing, as buoyant forces can inhibit vertical mixing. Likewise, 
the lack of buoyancy in our isothermal models may artificially enhance both the residual 
transport at the midplane and the coherent oscillations. In order to understand the effect of 
isothermality on these issues, we ran a single model at 32 zones per scale height including 
an ideal gas equation of state P = pkT / fi and an entropy equation, 

pT(— + u ■ Vs + ul — ) = r]{z)pof + Cp (V ■ uf + k^V^s, (6) 

where is a hyperdiffusivity of entropy, vPy = —3/2Qx is the Keplerian velocity profile, Q is 
the rotation rate of the shearing box, and all other symbols have their usual meanings. The 
entropy equation includes resistive and shock heating but not heating from the hyperdiffu- 
sivities (these terms are small, and do not significantly contribute to the thermal balance). 
We chose a ratio of specific heats 7 = 5/3. 

The boundary conditions on this run as on the others are periodic in all directions, which 
precludes the escape of heat. We choose this rather unrealistic setup to directly compare 
with our isothermal simulations, in order to demonstrate that our results are not due to the 
lack of buoyancy in the isothermal case. Because of the turbulent heating, the total energy 
in the box increases compared to the isothermal case. Although the turbulence is strongest 
in the surface layers, there is no temperature inversion in the vertical direction: the midplane 
remains hotter than the surface. The disk quasi-statically adjusts to new equilibrium density 
and pressure profiles as the energy increases. While these profiles are somewhat different 
from the standard Gaussian density profile expected for accretion disks with linear z-gravity 
profiles. Figure [5] shows that the main result of this difference is a factor ~ 2 increase in 
magnetic energy, and a factor ~ 1.8 increase in kinetic energy. 

The overall morphology of the flow is roughly similar to the isothermal model, with a 
reduced amplitude of motion in the dead zone but a similar large-scale structure. The stress 
profile across the dead zone is narrower but deeper in the ideal gas run than the isothermal 



- 11 - 



one because of the pressure confinement caused by turbulent heating in a periodic box absent 
coohng. We directly compare the Maxwell and Reynolds stresses for isothermal (labeled 
7 = 1) and ideal gas (7 = 5/3) in Figure 13 We find that inclusion of the buoyant forces 
reduces the amplitude of the vertical motions and increases their characteristic frequency 
to Wmax — 0.5f2o, still wcll below the first acoustic mode, as shown in Figure [14} However, 
the motions are not suppressed, even in this model in which turbulent heating is not at all 
balanced by cooling. A real disk in which cooling is allowed will fall between the limits of 
the isothermal model and this heating only model. Therefore, we argue that the dead zone 
motions represent a physical phenomenon even in the presence of buoyancy forces. 



4. Midplane Vortex Dynamics 



We see no vortices in our simulations. We must establish that this is a physical result, 
and not due to insufficient resolution. JG05 demonstrate that sufficient numerical resolution 
is necessary to sustain kinetic energy and angular momentum transport from vortices. In 
their 2D simulations, they find 128dx/H to be the critical value, raising the question of 
whether vortices can form in 2D at our resolutions {32dx/H or QAdx/H)! Furthermore, the 
typical radial size scale of a vortex in their simulations is roughly if, the entire width of our 
box. Can vortices form in such cramped quarters? 



4.1. Two-Dimensional Vortex Dynamics 



In order to clarify these issues, we ran four sets of 2D hydrodynamic simulations. First, 
we studied resolution. A comparison between Figure [2] and Figure 1 of JG05 shows that the 
sixth-order Pencil Code is significantly less diffusive than their second-order, Zeus-derived 
method, suggesting that we might be able to see sustained vortex activity at lower resolution. 
Therefore we ran a set of models using their domain with size = Ly = AH at resolutions 
of (16, 32, 64)(ix/if, all with Gaussian random velocity perturbations with the same initial 
perturbation amplitude, a = O.Sc^. Figure [TS] shows the kinetic energy and a-parameter, a = 
(Lu^Uy) for 2D. This Figure demonstrates clear convergence in these integrated quantities 
at 32dx/H. 

Next, we turn to the question of domain size. JG05 use a domain four times larger in 
radial extent than we do in our 3D simulations. Could this inhibit vortex formation in our 
MHD runs? To test this, we ran a set of 2D hydrodynamic models with the same domain 
as the midplane of our 3D simulations, Lx = IH and Ly = AH, again with resolutions of 
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(16, 32, 64)(ia;/if, seeded with perturbations of magnitude a = O.Scg. Figure 16 shows that 
vortex morphology is largely similar in both domains at three times, though there are more 
vortices present in the square domain. 



Figure 17 shows the resulting energy and a. Although the results are not as well 
converged as in the square domain, there is still evidence for vortex activity with lifetimes 
clearly greater than the MRI growth time tmri ~ i7 at a resolution of 32dx/B.. 

What, then, is the reason we do not see vortex activity in the 3D dead zone simulations? 
We have argued that dimensionality is not the limiting factor, as regions with lower Froude 
number that are more effectively 2D do not show more vorticity. It also appears clear that 
the 2D vortex activity found in previous simulations can be reproduced on domains like ours 
at the resolutions we use in our MHD runs, so resolution and grid size are also not limiting. 
Rather, the discriminant appears simply to be the strength of perturbations required to 
trigger long-lived vortices. We ran a third set of 2D hydrodynamic simulations, this time 
at a fixed resolution of 32o?a;//H, but with decreasing velocity perturbation magnitudes 
a = 0.8cs,0.5cs, and O.lc^. Figure [TS] shows that below a = 0.8, the kinetic energy and 
a drop precipitously after only a few orbits. (Note that we keep the domain size fixed at 
Lx = IH, Ly = 4H, unlike in the tests of JG05, in which the domain size was scaled linearly 
with the strength of the velocity dispersion.) We conclude that compressible vortices require 
strong a = O.SCg initial perturbations in order to survive for many orbits. 



Umurhan & Regev (2004) see perturbation energy sustained essentially indefinitely for 
Re = oo and decaying away only very slowly for Re = 50000, but their model is incompress- 
ible. They show this approximation to be rigorously valid for the small length scales they 
consider. However, interestingly, they find that the perturbation energy saturates at roughly 
< 10~^Eshear, whcrc Eghear IS the energy in the shear. They define the turbulent intensity of 
an incompressible shear flow as 

with u the disturbance velocity and Ush the shear velocity. We write Ush = qfix, integrate 
for the Keplerian case q = —3/2, and use our definition of (thermal) scale height H = 2cs/f2 
to arrive at 

with (Ma^)"'^ = Marms the RMS Mach number. Using their saturation value, this gives 
Ma^ms ~ 0.06{Lx/H), or about 0.06 over a domain of H. This is an order of magnitude 
below the levels necessary to incite vortices on the larger scales that we simulate in which 
compression becomes important. If we assume that small scale, non-linear, transient growth 



- 13 - 



instability scales with Reynolds number as found by Lesur & Longaretti (2005), not only is 
the incompressible, small-scale, transient growth irrelevant for angular momentum transport, 
its saturation level is so small that larger scale vortices in the compressible regime will not 
form at all. 



4.2. Long-lived axisymmetric structures 



Recently, Johansen et al. ( 2009 ) reported the existence of long-lived axisymmetric zonal 
flows in MRI turbulence, provided the radial width of the box is larger than ~ IH. These 
structures result from an inverse cascade of magnetic energy to the largest scales present in 
the box. This energy drives a large scale Maxwell stress, which in concert with the Coriolis 
force enforces the zonal flow. The end result are long-lived (~ 10s of orbits) pressure and 
density fluctuations. It is easy to imagine that these complex structures could have significant 
influence on motions in the dead zones, assuming there is sufficient magnetic energy to launch 
the zonal flows in the first place. Figure [T] implies that this criterion is unlikely to be met: 
even in the moderate sized dead zone models (e.g., Rcm = 30 and = 1); the total Maxwell 
stress is many orders of magnitude lower than that in the MRI active zones. However, given 



that Johansen et al. (2009) demonstrate that the zonal flows are subtle, stochastically driven. 



and non-linear, we find such a crude estimate to be unsatisfactory. 

Therefore, we re-ran our -Re a/ = 30 FS model with a radially-enlarged domain, running 
from —H < X < H . By averaging p over the entire y — z domain, including both active 
and dead zones, we recover the reported long-lived axisymmetric structures. Nevertheless, 



our results remain robust: there is no vortex formation, and figure 19 demonstrates that the 



midplane oscillations are quite similar between the larger domain and our fiducial case. 



5. Discussion 

We find a complete lack of coherent anticyclonic vortices in all of our stratified local 
simulations of dead zones, despite the presence of a vorticity source in the form of active 
zones. While our experiments do not tell us exactly what happens, it appears that even if 
the midplane were a pure 2D layer, large scale vortex formation would not occur because 
the amplitude of velocity perturbations driven in the dead zone by the active layers is well 



below that needed. Recent work by Lesur & Papaloizou (2009) suggests that vortices in 3D 



disks are always parametrically unstable, though the growth rates of such instabilities can 
be very slow. Even in the presence of an active vorticity source, though, we find no vortex 
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formation. 



Our simulations cover a smaller range of z than Barranco & Marcus (2005), and it is 



possible that the increased stratification present at higher z may change our results. To 



strengthen the results presented here, direct comparison with Barranco & Marcus (2005) 



including compressibility would be quite enlightening. Furthermore, almost all of our simu- 
lations use an isothermal relation for pressure and density. This forces density gradients and 
pressure gradients to coincide, eliminating the baroclinic term driving the vortex formation 



from breaking gravity waves in Barranco & Marcus (2005). Relaxing this requirement in 



our simulations may aid vortex formation at higher scale heights though our one adiabatic 
model without cooling did not show vortex formation either. 

However, we do find a large-scale, axisymmetric oscillation about the midplane in all 
dead zone models. This oscillation carries nearly all the kinetic energy of the dead zone. We 
are unable to associate it with a normal mode using a simple linear dispersion relation. A 
more detailed analysis including the effects of continuous stratification could shed consid- 
erable light on the situation. Regardless of such details, the low frequency (less than fig) 
suggests that the mode will not have much effect on dust sedimentation, which is typically 
dominated by random fluctuation with correlation times Tcorr — 0.15tor6, corresponding to 
frequencies ~ 7f2o (Fromang &; Papaloizou 2006). It could, in principle, cause a warp in the 



dust layer, which in turn could affect the gravitational stability (Goldreich & Ward 1973 



Johansen et al. 2007). Such phenomena can only be examined with a combined treatment 



of dust, self-gravity, and a dead zone, which we defer to future work. 

Finally, global dynamics may significantly alter vortex formation properties. A real 
protoplanetary disk is thin and quasi-2D for horizontal scales greater than the scale height. 
It is quite possible that robust, anticyclonic vortices may form within a disk, but our results 
suggest that a local mechanism is not responsible. The baroclinic instability (Klahr & 



Bodenheimer||2003| [Petersen et aL]|2007aP remains a candidate for forming such large scale 
vortices. 
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Fig. 1. — Maxwell stress as a function of height z for three different values of magnetic 
Reynolds number using the FS profile, and three different transition heights Zi^ using the 
tanh resistivity profile (Eq. [T]). At very high ri values, the stress is occasionally positive and 
therefore absent on this plot. 
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Fig. 2. — Time evolution of the velocity perturbation amplitude 5ux of an incompressible, 
shearing wave (upper panel) using physical viscosity for resolutions from 64^ to 128^. The 
analytic solution is given in black. Aliasing is present in the 64^ solutions but injects only 
trivial amounts of energy. The lower panel shows absolute value of the error, |e| = \5ux — 

^^analyticy ^^^^ ^^^^^ ^^^g 

error smoothly asymptotes to the analytic solution as the code 
hyperviscously damps the wave faster than aliasing can inject spurious energy. 
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Fig. 3. — Time evolution of the amplitude Suy of a compressible, shearing wave using 
hyperviscosity for resolutions from 32^ to 128^. The analytic solution is given in black. The 
lower panel shows the absolute value of the error, |e| = \6uy — ,5^i^"<^'y*«c|^ between the Pencil 
Code solution and a numerically integrated solution to the exact parabolic cylinder equation 
describing the wave. 
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Fig. 4. — Time evolution of the kinetic energy of a 3D, non-linear, incompressible, shearing 
wave using hyperviscosity for resolutions from 64'^ to 256^. A numerical integration of the 
exact wave differential equation is given in black and labeled "Analytic". The lower panel 
shows absolute value of error, |e| = \Ekin — E'^irf^^^'^l 
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Fig. 5. — From top to bottom: volume averaged vorticity (u), kinetic, and magnetic energy 
comparing isothermal and ideal gas Rcm = 30 runs. Isothermal runs at H = 64dx and 
H = 32dx are overplotted. The amount of vorticity increases with increasing resolution. 
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Fig. 7. — Internal Froude number versus height above the midplane for the FS runs. At 
Fr< 1, the flow becomes strongly stratified and effectively two-dimensional for vortices lying 
in the disk plane. 
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Fig. 8. — Vorticity in horizontal {x-y) planes as a function of height above and below the 
midplane for the Rcm = 3 run. The effectively 2D nature of the flow does not appear to 
create regions in which coherent vortices can form. 
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Fig. 9. — Mean square velocities in each direction for the FS runs. In all cases with a 
dead zone, the upper three curves show the active zone and the lower three the dead zone 
velocities. Note the transition in dominant component from y velocity to z velocity from 
active to dead in all cases. 
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Fig. 10. — Slices through the midplane of density (color) with arrows following the in-plane 
{x-y) velocity. Each density image is scaled to its own minimum and maximum to emphasize 
morphology. 
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Fig. 11. — Slices through the y = (radial- vertical) plane of density (color) with arrows 
following the in-plane {x-z) velocity. Each density image is scaled to its own minimum and 
maximum to emphasize morphology. Note the vertical circulation in dead zones viewed in 
this plane. 
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Fig. 12. — Temporal power spectra for Ux, Uy, p. Acoustic modes are present in all runs at 
frequencies > IQ. Dead zone models also show a low-frequency inertial oscillation apparently 
driven by the active zones. 
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Fig. 14. — Temporal power spectra of u^, Uy, p for Rcm = 30 with ideal gas equation of 
state (labeled 7 = 5/3) and isothermal (7 = 1) runs overplotted. Both inertial and acoustic 
waves are significantly damped but still present in the 7 = 5/3 case. The large amplitude 
mode at vuQ ~ 0.2 in the isothermal case remains present, though at reduced amplitude and 
increased wavenumber vuQ ~ 0.5. The normalization difference between the two spectra is 
due to the extra heat retained in the non-isothermal case. 
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Fig. 15. — Kinetic energy and a for a 2D domain with = Ly = AH at three different 
resolutions, (16, 32, 64)(ix/if, using initial velocity perturbations of magnitude a = O.Scg. 
Energy and a are clearly well converged at a resolution of 32dx/H. 
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Fig. 16. — Potential vorticity of the entire x-y plane at times tQ = 15.7,31.4,47.1 (left to 
right) upper panel A AH square domain. Blue tones indicate negative potential vorticity, red 



positive. Vortex formation is similar to that found in Johnson & Gammie (2005b). lower 



panel A domain IH x AH. In both panels, colors are scaled to the minimum and maximum 
of each image to highlight vortex morphology. 
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Fig. 17. — Kinetic energy and a for a domain with = IH, and Ly = AH, at three different 
resolutions, (16, 32, 64)(ix/if, using initial velocity perturbations of magnitude cr = O.Scs- 
Energy and a are not as well converged as in Figure 15, but sustained vortex activity does 
seem to occur even at 32 zones/H. 




Fig. 18. — Kinetic energy and a for a domain with = IH and Ly = AH at resolution 
of 32dx/H, using three different initial velocity dispersions, a = (0.1, 0.5, 0.8)cs. Sustained 
vortex activity clearly requires strong initial perturbations. 
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Fig. 19. — Temporal power spectra of u^, Uy, p for Rcm = 30 with the standard box width 
{~0.5H < X < 0.5H) and a box twice as wide {—H < x < H). The wider box features the 



large axisymmetric pressure bumps reported by Johansen et al.| (2009), but this has little 
effect on the oscillations reported in this paper. 
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Table 1. Directional normalized enstrophy for dead and active zones. 







Active 






Dead 




Run 














64Rinf 


0.31 


0.34 


0.35 








64R100 


0.32 


0.34 


0.35 


0.11 


0.69 


0.20 


64R30 


0.35 


0.36 


0.29 


0.07 


0.79 


0.14 


64R3 


0.31 


0.41 


0.28 


0.06 


0.68 


0.26 


32z0.5 


0.26 


0.48 


0.27 


0.10 


0.68 


0.22 


32zl.O 


0.27 


0.46 


0.26 


0.07 


0.76 


0.17 



